function [outSt,logL]=kfilterRegSplit(parvec,funcmod,Y,trainvec,solveopt,addsol)
% ====================================================================
% ====================================================================
%% 1. Model solution
% Matrices of second sample stored in structure second 
[G,R,C,eu,SDX,Z,structOne,~,structTwo]=feval(funcmod,parvec,solveopt,addsol);
% A. Non-existence in either sample 
if isequal(eu,[1;1])==0 
    return
end
[T,ny]=size(Y);
%% 2. Demeaning using the split sample 
%
%  First Sample (C1)
if any(C~=0)==1
    Y(1:addsol.rowBegSecond-1,:)=...
        Y(1:addsol.rowBegSecond-1,:)-repmat((Z*C)',[addsol.rowBegSecond-1 1]);
end 
% Second sample (second.second.C)
if any(structOne.C~=0)==1;
    Y(addsol.rowBegSecond:end,:)=...
        Y(addsol.rowBegSecond:end,:)-repmat((structOne.Z*structOne.C)',[(T-addsol.rowBegSecond+1) 1]);
end 
Y=Y';
%% 3. Forward filter 

%% 3.1. Dimensions and storage 
 
ns=size(G,1); 
nx=size(SDX,1); 

vt=zeros(ny,T);
finvt=zeros(ny,ny,T);
kpartg=zeros(ns,ny,T);
logLnc=zeros(T,1); 
yfor =zeros(ny,T);

% Matrices with one additional entry (initialization)
% to recover observables 
stt=zeros(ns,T+1);
ptt=zeros(ns,ns,T+1);

%% 3.2 Initialization 

stt(:,1)=zeros(ns,1);
ptt(:,:,1)=feval(@lyapunov_symm,G,R*(SDX')*(R*(SDX')'));

%% 3.3 Start Forward Filter using KF_DK 
% First Part
for ii=1:addsol.rowBegSecond-1;
    yfor(:,ii) =Z*(G*stt(:,ii));
    [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vt(:,ii),finvt(:,:,ii),kpartg(:,:,ii)]...
        =feval(@kf_dk,Y(:,ii),Z,stt(:,ii),ptt(:,:,ii),G,R*(SDX'));
end
% Second part
for ii=addsol.rowBegSecond:Nobs;
    yfor(:,ii) =structOne.Z*(structOne.G*stt(:,ii));
    [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vt(:,ii),finvt(:,:,ii),kpartg(:,:,ii)]=...
        feval(@kf_dk,Y(:,ii),structOne.Z,stt(:,ii),ptt(:,:,ii),structOne.G,...
        (structOne.R)*(structOne.SDX'));
end

%% 4. Likelihood with Integration Constant 
 
logLnc=logLnc(trainvec(1):trainvec(2)); 
logL=-0.5*log(2*pi)*(ny*T)+sum(logLnc);

%% 5. Truncate filters and obtain initial observations 
yferr=vt';
yfor =yfor'; 
stt=stt(:,2:end); 
ptt=ptt(:,:,2:end); 

%% 6. Disturbance smoother with TV matrices 
% Obtain the Innovations using a disturbance smoother 
etamat=zeros(NX,Nobs); 
smooth_st=zeros(NS,Nobs); 
rmat=zeros(NS,Nobs); 

%% 6.1 Initialize RSTAR & start at t=Nobs 
rstar=zeros(NS,1); 
[rstar,etamat(:,end)]=feval(@smoothdis,rstar,(structOne.SDX')*structOne.SDX,...
    structOne.R',structOne.Z',finvt(:,:,end),zeros(ns),vt(:,end));
smooth_st(:,end)=stt(:,end)+ptt(:,:,end)*rstar; 
rmat(:,end)=rstar; 

%% 6.2 Begin Backward recursion 
Ztr=structOne.Z';
Qf = (structOne.SDX')*structOne.SDX;
Rf = structOne.R;
Gf = structOne.G;

for ii=(Nobs-1):-1:1;
    % Recall addsol.rowBegSecond is the beginning of the Second Sample 
    % This break must occur exactly when second sample begins
    if ii==addsol.rowBegSecond-1;
        Ztr=Z';
        Qf=(SDX'*SDX);
        Rf=R; 
    end
    [rstar,etamat(:,ii)]=feval(@smoothdis,rstar,Qf,(Rf'),...
        Ztr,finvt(:,:,ii),((Gf-Gf*kpartg(:,:,ii)*Zf)'),vt(:,ii));
    smooth_st(:,ii)=stt(:,ii)+ptt(:,:,ii)*rstar;
    rmat(:,ii)=rstar; 
    % Recall addsol.rowBegSecond -1 is the end of the First Sample
    % This break must occur exactly when first sample ends
    if ii==addsol.rowBegSecond-1; 
        Gf=G;          
    end         
end
a0 = P0*(Gf')*rstar;   % Note: this is only correct in the case where a0 = zeros(ns,1) 
etamat=(etamat)';
smooth_st=(smooth_st)'; 
stt  =stt'; 


%% 7. Check Smoother 
% Check that Smooth States are identical if using disturbance smoother (above) vs. state smoother (below) 
% and if can also recover the observables
outcount.aold=azero;
outcount.Pold=Pzero;
tol=1e-5; 
% State at time zero give by GG1*azero+Pzero*rstar;
[yCheck,smoothCheck]=feval(@kfilterRegSplitSimulation,a0,etamat'); 
disp('Max Discrepancy Smooth vs. Actual Data'); 
maxdifY=feval(@comparemat,yCheck,Y'); 
disp('Max Discrepancy State and Innovation Smoother'); 
maxdifS=feval(@comparemat,smoothCheck,smooth_st); 
if maxdifY > tol || maxdifS > tol 
    error('Smoother discrepancy exceeds tolerance') 
end 

%% 8. Initial Variance and State per shock (used below for the counterfactual decompositions)
% vInitialPershock: Initial Variance, decomposed per shock 
% sOnePerShock, State at time 1 decomposed, decomposed per shock
vInitialPerShock=zeros(ns,ns,nx);
sZeroSmoothPerShock=zeros(ns,nx); 
sOneSmoothPerShock =zeros(ns,nx); 
for ii=1:nx
    vInitialPerShock(:,:,ii)=feval(@lyapunov_symm,G,R(:,ii)*SDX(ii,ii)*SDX(ii,ii)*(R(:,ii)'));
    sZeroSmoothPerShock(:,ii) = vInitialPerShock(:,:,ii)*(G')*rstar;
    etatemp = zeros(nx,1);
    etatemp(ii) = etamat(1,ii);
    sOneSmoothPerShock(:,ii)=G*sZeroSmoothPerShock(:,ii)...
        +R*etatemp;
end
maxdifSzero=comparemat(sum(sZeroSmoothPerShock,2),a0); 
maxdifSone =comparemat(sum(sOneSmoothPerShock,2),smooth_st(1,:)' ); 
if maxdifSzero > tol; error('Cannot decompose sZeroPerShock'); end 
if maxdifSone > tol; error('Cannot decompose sOnePerShock'); end
%% 9. Counterfactual Decomposition 
% Generate states by feeding each shock at a time, ensure that it recovers
% the original state 
disp('Begin Counterfactual')
countStates=zeros(T,ns,nx);
countObs   =zeros(T,ny,nx);
for ii=1:nx
    etaTemp=zeros(T,nx); 
    etaTemp(:,ii)=etamat(:,ii); 
    [countObs(:,:,ii),countStates(:,:,ii)]=feval(@kfilterRegSplitSimulation,sZeroSmoothPerShock(:,ii),etaTemp');
end 
disp('Maximum Discrepancy Counterfactual States and Smooth States') 
maxdifCount=comparemat(sum(countStates,3),smooth_st); 
if maxdifCount > tol; error('Counterfactual States do not recover smooth states'); end 

outSt.filteredStates=stt; 
outSt.smoothStates=smooth_st; 
outSt.innovations=etamat; 
outSt.countStates=countStates; 
outSt.countObs=countObs; 
outSt.decompInitialState=sOneSmoothPerShock; 
outSt.forecastObs=yfer; 
outSt.forecastError=yferr; 

   
%% Subroutine kfilterRegSplitSimulation allows to simulate the model  
% Inputs 
% sInitial: [NS 1] Initial State 
% innovMat: [Nx Nobx] matrix of innovations 
% By being a sub-routine it has access to all variables defined above 
% Be-careful not to use an index (ii,jj) used above or to repeat variable
% names
    function [ySim,sSim]=kfilterRegSplitSimulation(sInitial,innovMat) 
        if ~isequal(size(innovMat),[nx T]) 
            error('Input innovMat must be [NX Nobs]') 
        end 
        sSim=zeros(ns,T);
        ySim=zeros(size(Y));
        sSim(:,1)=G*sInitial+R*innovMat(:,1);
        ySim(:,1)=Z*sSim(:,1);
        for kk=2:addsol.rowBegSecond-1;
            sSim(:,kk)=G*sSim(:,kk-1)+R*innovMat(:,kk);
            ySim(:,kk)=Z*sSim(:,kk);
        end
        for kk=addsol.rowBegSecond:T;
            sSim(:,kk)=structOne.G*sSim(:,kk-1)+structOne.R*(innovMat(:,kk));
            ySim(:,kk)=structOne.Z*sSim(:,kk);
        end       
        ySim=ySim'; 
        sSim=sSim'; 
    end
%% End of File
end 